---
title: "Research and Politics Replication Code: The Treatment Works"
author: "Bryan Wilcox-Archuleta"
date: "Feb 4, 2018"
output: pdf_document
---

This section loads required libraries. 

```{r, message=FALSE, warning=FALSE}
library(tidyverse)
library(gridExtra)
library(simcf) # must be installed from http://faculty.washington.edu/cadolph/index.php?page=60
library(stargazer)
library(RColorBrewer)
library(MASS)
col <- brewer.pal(3,"Dark2")
ggplot.text <- element_text(size = 12)

```

# Read in the data

```{r, eval=FALSE}
data <- read_csv("replication_data.csv") # eval set to false
```


# Loess by group 

This makes Figures 1-4 that are in the manuscript. 

```{r, message=FALSE, warning=FALSE}
# all groups 
latinoLoess <- ggplot(data, aes(y=latinoHood, x=latinoTract)) + 
  geom_point(pch=1, col="gray35") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Latino\n Neighborhood",y="Perceived composition", x="Proportion Latino") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),
                     labels=c("Not At\n All", "Mostly", "Almost\n Entirely"))

blackLoess <- ggplot(data, aes(y=blackHood, x=blackTract)) + 
  geom_point(pch=1, col="gray35") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Black\n Neighborhood", y="",x="Proportion Black") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

whiteLoess <- ggplot(data, aes(y=whiteHood, x=whiteTract)) + 
  geom_point(pch=1, col="gray35") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="White\n Neighborhood", y="",x="Proportion White") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

asianLoess <- ggplot(data, aes(y=asianHood, x=asianTract)) + 
  geom_point(pch=1, col="gray35") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Asian\n Neighborhood", y="",x="Proportion Asian") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

mixedLoess <- ggplot(data, aes(y=mixedHood, x=escore_tract)) + 
  geom_point(pch=1, col="gray35") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Mixed\n Neighborhood", y="",x="Entropy") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 


# Figure 1
grob <- gridExtra::arrangeGrob(latinoLoess, blackLoess, whiteLoess, asianLoess, mixedLoess,nrow=1, 
                               widths = c(.8,.7,.7,.7,.7))
ggsave(filename ="loess_all_groups.png",plot = grob, width = 12, height = 6)


# Figure 2
# for Latinos
latinoData <- data %>% filter(latino == 1)
latinoLoess <- ggplot(latinoData, aes(y=latinoHood, x=latinoTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="lm", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Latino\n Neighborhood",y="Perceived composition", x="Proportion Latino") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),
                     labels=c("Not At\n All", "Mostly", "Almost\n Entirely"))


blackLoess <- ggplot(latinoData, aes(y=blackHood, x=blackTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Black\n Neighborhood", y="",x="Proportion Black") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", ""))

whiteLoess <- ggplot(latinoData, aes(y=whiteHood, x=whiteTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="White\n Neighborhood", y="",x="Proportion White") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 


asianLoess <- ggplot(latinoData, aes(y=asianHood, x=asianTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Asian\n Neighborhood", y="",x="Proportion Asian") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

mixedLoess <- ggplot(latinoData, aes(y=mixedHood, x=escore_tract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Mixed\n Neighborhood", y="",x="Entropy") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

grob <- gridExtra::arrangeGrob(latinoLoess, blackLoess, whiteLoess, 
                               asianLoess,mixedLoess,nrow=1, widths = c(.8,.7,.7,.7,.7))
ggsave(filename ="loess_latino.png",plot = grob, width = 12, height = 6)

# Figure 3
# for Blacks
blackData <- data %>% filter(black == 1)
latinoLoess <- ggplot(blackData, aes(y=latinoHood, x=latinoTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Latino\n Neighborhood",y="Perceived composition", x="Proportion Latino") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),
                     labels=c("Not At\n All", "Mostly", "Almost\n Entirely"))

blackLoess <- ggplot(blackData, aes(y=blackHood, x=blackTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Black\n Neighborhood", y="",x="Proportion Black") +
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", ""))

whiteLoess <- ggplot(blackData, aes(y=whiteHood, x=whiteTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="White\n Neighborhood", y="",x="Proportion White") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

asianLoess <- ggplot(blackData, aes(y=asianHood, x=asianTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Asian\n Neighborhood", y="",x="Proportion Asian") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

mixedLoess <- ggplot(blackData, aes(y=mixedHood, x=escore_tract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Mixed\n Neighborhood", y="",x="Entropy") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 


grob <- gridExtra::arrangeGrob(latinoLoess, blackLoess, whiteLoess, 
                               asianLoess,mixedLoess,nrow=1, 
                               widths = c(.8,.7,.7,.7,.7))

ggsave(filename ="loess_black.png",plot = grob, width = 12, height = 6)

# Figure 4
# for whites
whiteData <- data %>% filter(white == 1)
latinoLoess <- ggplot(whiteData, aes(y=latinoHood, x=latinoTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Latino\n Neighborhood",y="Perceived composition", x="Proportion Latino") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),
                     labels=c("Not At\n All", "Mostly", "Almost\n Entirely"))

blackLoess <- ggplot(whiteData, aes(y=blackHood, x=blackTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="Black\n Neighborhood", y="",x="Proportion Black") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", ""))

whiteLoess <- ggplot(whiteData, aes(y=whiteHood, x=whiteTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + 
  labs(title="White\n Neighborhood", y="",x="Proportion White") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

asianLoess <- ggplot(whiteData, aes(y=asianHood, x=asianTract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Asian\n Neighborhood", y="",x="Proportion Asian") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

mixedLoess <- ggplot(whiteData, aes(y=mixedHood, x=escore_tract)) + 
  geom_point(pch=1, col="gray25") + 
  geom_smooth(method="loess", lwd=1.5, col = "black") + theme_bw() + theme_bw() + 
  labs(title="Mixed\n Neighborhood", y="",x="Entropy") + 
  scale_y_continuous(limits=c(0,2),breaks=c(0,1,2),labels=c("", "", "")) 

grob <- gridExtra::arrangeGrob(latinoLoess, blackLoess, 
                               whiteLoess, asianLoess,mixedLoess,nrow=1, 
                               widths = c(.8,.7,.7,.7,.7))
ggsave(filename ="loess_white.png",plot = grob, width = 12, height = 6)
```

# Model

This code makes model 1. 

```{r, message=FALSE, warning=FALSE}
# Table 1 
latino_model <- as.factor(latinoHood) ~ latinoTract + edu + income + age + female + 
  latino + black + democrat + own_home + pocketbookEval + cmps2012
black_model <- as.factor(blackHood) ~ blackTract + edu + income + age + female +
  latino + black + democrat + own_home + pocketbookEval + cmps2012
white_model <- as.factor(whiteHood) ~ whiteTract + edu + income + age + female +
  latino + black + democrat + own_home + pocketbookEval + cmps2012
asian_model <- as.factor(asianHood) ~ asianTract + edu + income + age + female +
  latino + black + democrat + own_home + pocketbookEval + cmps2012
mixed_model <- as.factor(mixedHood) ~ escore_tract + edu + income + age + female +
  latino + black + democrat + own_home + pocketbookEval + cmps2012


latinoResult <- polr(latino_model, data = data, method = "probit", Hess = TRUE)
blackResult <- polr(black_model, data = data, method = "probit", Hess = TRUE)
whiteResult <- polr(white_model, data = data, method = "probit", Hess = TRUE)
asianResult <- polr(asian_model, data = data, method = "probit", Hess = TRUE)
mixedResult <- polr(mixed_model, data = data, method = "probit", Hess = TRUE)

labs <- c("\\% Latino", "\\%Black", "\\% White","\\% Asian American","Education",
          "Income","Age","Female","Latino","Black","Democrat","Home Owner", 
          "Pocketbook Evaluation","2012 CMPS")

stargazer(latinoResult, blackResult, whiteResult, asianResult,
          title = "\\textbf{Relationship between Census Measured and Perceived Composition}",
          label = "table1",
          covariate.labels = labs,
          star.cutoffs = c(0.05, 0.01, 0.001), 
          dep.var.labels = c('Latino', 'Black', 'White', 'Asian American'),
          model.numbers          = FALSE,
          out = "table_1.tex")
```

# Appendix

This code makes the figures for the appendix. 

```{r, message=FALSE, warning=FALSE}

latinoData <- data %>% dplyr::select(latinoHood, latinoTract, blackHood, blackTract, 
                                     asianHood, asianTract, whiteHood, whiteTract, 
                                     mixedHood, escore_tract, cmps2012, res_length, 
                                     home_own)

#model 
latinoTractModel <- as.factor(latinoHood) ~ latinoTract + home_own + res_length + cmps2012
latinoTractData <- extractdata(latinoTractModel, data=latinoData, na.rm = TRUE)
latinoTractResult <- polr(latinoTractModel, data = latinoTractData,
                          method = "probit", Hess = TRUE)


blackTractModel <- as.factor(blackHood) ~ blackTract + home_own + res_length + cmps2012
blackTractData <- extractdata(blackTractModel, data=latinoData, na.rm = TRUE)
blackTractResult <- polr(blackTractModel, data = blackTractData, 
                         method = "probit", Hess = TRUE)


asianTractModel <- as.factor(asianHood) ~ asianTract + home_own + res_length + cmps2012
asianTractData <- extractdata(asianTractModel, data=latinoData, na.rm = TRUE)
asianTractResult <- polr(asianTractModel, data = asianTractData, 
                         method = "probit", Hess = TRUE)


whiteTractModel <- as.factor(whiteHood) ~ whiteTract + home_own + res_length + cmps2012
whiteTractData <- extractdata(whiteTractModel, data=latinoData, na.rm = TRUE)
whiteTractResult <- polr(whiteTractModel, data = whiteTractData, method = "probit", Hess = TRUE)


mixedTractModel <- as.factor(mixedHood) ~ escore_tract + home_own + res_length + cmps2012
mixedTractData <- extractdata(mixedTractModel, data=latinoData, na.rm = TRUE)
mixedTractResult <- polr(mixedTractModel, data = mixedTractData, 
                         method = "probit", Hess = TRUE)


############################
latinoTractSims <- mvrnorm(1000,c(latinoTractResult$coefficients, 
                                  latinoTractResult$zeta), 
                           vcov(latinoTractResult))

latinoTractModelSim <- latinoHood ~ latinoTract + home_own + res_length + cmps2012

latinoTractSeq <- seq(0,1,.1)
latinoTractMatrix <- cfMake(latinoTractModelSim, data = latinoTractData, 
                            nscen = length(latinoTractSeq))

for(i in 1:length(latinoTractSeq)){
  latinoTractMatrix <- cfChange(latinoTractMatrix, "latinoTract", 
                                x=latinoTractSeq[i], scen = i)
}
latinoTractResultEV <- oprobitsimev(latinoTractMatrix, latinoTractSims, 
                                    constant = NA, cat = 3, ci=0.95)

###########################
# Black 

blackTractSims <- mvrnorm(1000,c(blackTractResult$coefficients, 
                                 blackTractResult$zeta), 
                          vcov(blackTractResult))

blackTractModelSim <- blackHood ~ blackTract + home_own + res_length + cmps2012

blackTractSeq <- seq(0,1,.1)
blackTractMatrix <- cfMake(blackTractModelSim, data = blackTractData, 
                           nscen = length(blackTractSeq))

for(i in 1:length(blackTractSeq)){
  blackTractMatrix <- cfChange(blackTractMatrix, "blackTract", 
                               x=blackTractSeq[i], scen = i)
}
blackTractResultEV <- oprobitsimev(blackTractMatrix, 
                                   blackTractSims,constant = NA, cat = 3, ci=0.95)


#################################
# Asian

asianTractSims <- mvrnorm(1000,c(asianTractResult$coefficients, 
                                 asianTractResult$zeta), 
                          vcov(asianTractResult))
asianTractModelSim <- asianHood ~ asianTract + home_own + res_length + cmps2012

asianTractSeq <- seq(0,1,.1)
asianTractMatrix <- cfMake(asianTractModelSim, data = asianTractData,
                           nscen = length(asianTractSeq))

for(i in 1:length(asianTractSeq)){
  asianTractMatrix <- cfChange(asianTractMatrix, "asianTract", 
                               x=asianTractSeq[i], scen = i)
}
asianTractResultEV <- oprobitsimev(asianTractMatrix, asianTractSims,
                                   constant = NA, cat = 3, ci=0.95)

############
# White

whiteTractSims <- mvrnorm(1000,c(whiteTractResult$coefficients, 
                                 whiteTractResult$zeta), 
                          vcov(whiteTractResult))

whiteTractModelSim <- whiteHood ~ whiteTract + home_own + res_length + cmps2012

whiteTractSeq <- seq(0,1,.1)
whiteTractMatrix <- cfMake(whiteTractModelSim, data = whiteTractData, 
                           nscen = length(whiteTractSeq))

for(i in 1:length(whiteTractSeq)){
  whiteTractMatrix <- cfChange(whiteTractMatrix, "whiteTract", 
                               x=whiteTractSeq[i], scen = i)
}
whiteTractResultEV <- oprobitsimev(whiteTractMatrix, whiteTractSims,
                                   constant = NA, cat = 3, ci=0.95)

##############
# Mixed

mixedTractSims <- mvrnorm(1000,c(mixedTractResult$coefficients, 
                                 mixedTractResult$zeta), 
                          vcov(mixedTractResult))

mixedTractModelSim <- mixedHood ~ escore_tract + home_own + res_length + cmps2012

mixedTractSeq <- seq(min(latinoData$escore_tract,na.rm = T),
                     max(latinoData$escore_tract,na.rm = T),.1)
mixedTractMatrix <- cfMake(mixedTractModelSim, 
                           data = latinoData, nscen = length(mixedTractSeq))

for(i in 1:length(mixedTractSeq)){
  mixedTractMatrix <- cfChange(mixedTractMatrix, "escore_tract", 
                               x=mixedTractSeq[i], scen = i)
}
mixedTractResultEV <- oprobitsimev(mixedTractMatrix, mixedTractSims, 
                                   constant = NA, cat = 3, ci=0.95)

mhyp <- seq(1, length(mixedTractSeq), 1)
mixedTractDF <- data.frame(x = mhyp, 
                           y = mixedTractResultEV$pe, 
                           ymin = mixedTractResultEV$lower,
                           ymax = mixedTractResultEV$upper)


xhyp <- seq(1, length(whiteTractSeq), 1)
latinoResultsDF<- data.frame(x = xhyp, 
                           latino.y.white = whiteTractResultEV$pe, 
                           latino.ymin.white = whiteTractResultEV$lower,
                           latino.ymax.white = whiteTractResultEV$upper,
                           latino.y.asian = asianTractResultEV$pe, 
                           latino.ymin.asian = asianTractResultEV$lower,
                           latino.ymax.asian = asianTractResultEV$upper,
                           latino.y.black = blackTractResultEV$pe, 
                           latino.ymin.black = blackTractResultEV$lower,
                           latino.ymax.black = blackTractResultEV$upper,
                           latino.y.latino = latinoTractResultEV$pe, 
                           latino.ymin.latino = latinoTractResultEV$lower,
                           latino.ymax.latino = latinoTractResultEV$upper)

xhyp <- seq(1, length(whiteTractSeq), 1)
latinoPlotLatino <- ggplot(latinoResultsDF, aes(x = xhyp, y = latino.y.latino.1)) + 
  geom_line(color = col[1]) + 
  geom_ribbon(aes(ymin=latino.ymin.latino.1,ymax=latino.ymax.latino.1),
              alpha=0.3, fill = col[1]) + 
  geom_line(aes(x=xhyp, y=latino.y.latino.2), color = col[2]) + 
  geom_line(aes(x=xhyp, y=latino.y.latino.3), color = col[3]) + 
  geom_ribbon(aes(ymin=latino.ymin.latino.2,ymax=latino.ymax.latino.2),
              alpha=0.3, fill = col[2]) +
  geom_ribbon(aes(ymin=latino.ymin.latino.3,ymax=latino.ymax.latino.3),
              alpha=0.3, fill = col[3]) +
  theme_bw() + 
  labs(y="Predicted Probability", x="% Latino in Tract", title="Latino Neighborhood") + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_discrete(breaks = c(1,3,5,7,9,11), 
                   labels = c("0%","20%","40%","60%","80%","100%")) + 
  annotate("text", label = "Non-Latino", x = 4, y = 0.8, 
           color = col[1], size=4) + 
  annotate("text", label = "Mostly Latino", x = 3, y = 0.15, 
           color = col[2], size=4) + 
  annotate("text", label = "Almost Entirely \n Latino", x = 9.5, y = 0.6, 
           color = col[3], size=4) + 
  theme(text = ggplot.text)

latinoPlotBlack <- ggplot(latinoResultsDF, aes(x = xhyp, y = latino.y.black.1)) + 
  geom_line(color = col[1]) + 
  geom_ribbon(aes(ymin=latino.ymin.black.1,ymax=latino.ymax.black.1),
              alpha=0.3, fill = col[1]) + 
  geom_line(aes(x=xhyp, y=latino.y.black.2), color = col[2]) + 
  geom_line(aes(x=xhyp, y=latino.y.black.3), color = col[3]) + 
  geom_ribbon(aes(ymin=latino.ymin.black.2,ymax=latino.ymax.black.2),
              alpha=0.3, fill = col[2]) +
  geom_ribbon(aes(ymin=latino.ymin.black.3,ymax=latino.ymax.black.3),
              alpha=0.3, fill = col[3]) +
  theme_bw() + labs(y="Predicted Probability", x="% Black in Tract", 
                    title="Black Neighborhood") + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_discrete(breaks = c(1,3,5,7,9,11), 
                   labels = c("0%","20%","40%","60%","80%","100%")) + 
  theme(axis.text.x = ggplot.text) + 
  annotate("text", label = "Non-Black", x = 4, y = 0.8, 
           color = col[1], size=4) + 
  annotate("text", label = "Mostly Black", x = 3, y = 0.05, 
           color = col[2], size=4) + 
  annotate("text", label = "Almost Entirely \n Black", x = 8.5, y = 0.4, 
           color = col[3], size=4) + 
  theme(text = ggplot.text)

latinoPlotWhite <- ggplot(latinoResultsDF, aes(x = xhyp, y = latino.y.white.1)) + 
  geom_line(color = col[1]) + 
  geom_ribbon(aes(ymin=latino.ymin.white.1,ymax=latino.ymax.white.1),
              alpha=0.3, fill = col[1]) + 
  geom_line(aes(x=xhyp, y=latino.y.white.2), color = col[2]) + 
  geom_line(aes(x=xhyp, y=latino.y.white.3), color = col[3]) + 
  geom_ribbon(aes(ymin=latino.ymin.white.2,ymax=latino.ymax.white.2),
              alpha=0.3, fill = col[2]) +
  geom_ribbon(aes(ymin=latino.ymin.white.3,ymax=latino.ymax.white.3),
              alpha=0.3, fill = col[3]) +
  theme_bw() + labs(y="Predicted Probability", x="% White in Tract",
                    title="White Neighborhood") + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_discrete(breaks = c(1,3,5,7,9,11), 
                   labels = c("0%","20%","40%","60%","80%","100%")) + 
  theme(axis.text.x = ggplot.text) + 
  annotate("text", label = "Non-White", x = 3.5, y = 0.8, 
           color = col[1], size=4) + 
  annotate("text", label = "Mostly White", x = 3, y = 0.15, 
           color = col[2], size=4) + 
  annotate("text", label = "Almost Entirely \n White", x = 9.5, y = 0.55, 
           color = col[3], size=4) + 
  theme(text = ggplot.text)

latinoPlotAsian <- ggplot(latinoResultsDF, aes(x = xhyp, y = latino.y.asian.1)) + 
  geom_line(color = col[1]) + 
  geom_ribbon(aes(ymin=latino.ymin.asian.1,ymax=latino.ymax.asian.1),
              alpha=0.3, fill = col[1]) + 
  geom_line(aes(x=xhyp, y=latino.y.asian.2), color = col[2]) + 
  geom_line(aes(x=xhyp, y=latino.y.asian.3), color = col[3]) + 
  geom_ribbon(aes(ymin=latino.ymin.asian.2,ymax=latino.ymax.asian.2),
              alpha=0.3, fill = col[2]) +
  geom_ribbon(aes(ymin=latino.ymin.asian.3,ymax=latino.ymax.asian.3),
              alpha=0.3, fill = col[3]) +
  theme_bw() + labs(y="Predicted Probability", x="% Asian-American in Tract", 
                    title="Asian American Neighborhood") + 
  scale_y_continuous(limits=c(0,1)) + 
  scale_x_discrete(breaks = c(1,3,5,7,9,11), 
                   labels = c("0%","20%","40%","60%","80%","100%")) + 
  theme(axis.text.x = ggplot.text) + 
  annotate("text", label = "Non-Asian\n American", x = 3.0, y = 0.8, 
           color = col[1], size=4) + 
  annotate("text", label = "Mostly \n Asian American", x = 3, y = 0.25, 
           color = col[2], size=4) + 
  annotate("text", label = "Almost Entirely \n Asian American", x = 9, y = 0.85, 
           color = col[3], size=4)+ 
  theme(text = ggplot.text)

# this is figure SI 0.5 
pdf(file="mpsa_appendix1.pdf",width = 12,height = 12)
grid.arrange(latinoPlotLatino,latinoPlotBlack,latinoPlotWhite,latinoPlotAsian,nrow=2)
dev.off()

mhyp <- seq(1, length(mixedTractSeq), 1)

# this is figure SI 0.6
PlotMixed <- ggplot(mixedTractDF, aes(x=mhyp,y=y.1)) + 
  geom_line(color = col[1])+ 
  geom_ribbon(aes(x=mhyp,ymin=ymin.1,ymax=ymax.1),
              alpha=0.3, fill = col[1]) + 
  geom_line(aes(x=mhyp, y=y.2), color = col[2]) + 
  geom_line(aes(x=mhyp, y=y.3), color = col[3]) + 
  geom_ribbon(aes(x=mhyp,ymin=ymin.2,ymax=ymax.2),alpha=0.3, 
              fill = col[2]) +
  geom_ribbon(aes(x=mhyp,ymin=ymin.3,ymax=ymax.3),alpha=0.3, 
              fill = col[3]) +
  theme_bw() + labs(y="Predicted Probability", x="Entropy in Tract",
                    title="Mixed Neighborhood") + 
  scale_y_continuous(limits=c(-0,1)) +
  scale_x_continuous(limits=c(1,13), breaks = mhyp, 
                     labels = round(mixedTractSeq,2)) + 
  annotate("text", label = "Non-mixed", x = 4.5, y = 0.8, 
           color = col[1], size=4) + 
  annotate("text", label = "Mostly \n mixed", x = 2, y = 0.25,
           color = col[2], size=4) + 
  annotate("text", label = "Almost Entirely \n mixed", x = 10, y = 0.08, 
           color = col[3], size=4) + 
  theme(text = element_text(size = 12))

pdf(file="mixed_pooled.pdf",width = 10,height = 10)
PlotMixed
dev.off()
```

# Residential Length Interaction

```{r, message=FALSE, warning=FALSE}
modelTract_latino <- as.factor(latinoHood) ~ latinoTract*res_length + 
  edu + income + age + female + latino + black + democrat + own_home + pocketbookEval
tractData_latino <- extractdata(modelTract_latino, data = data, na.rm = TRUE)
resTract_latino <- polr(modelTract_latino, data = tractData_latino, 
                        method = "probit", Hess = TRUE)
summary(resTract_latino)

modelTract_black <- as.factor(blackHood) ~ blackTract*res_length + 
  edu + income + age + female +latino + black + democrat + own_home + pocketbookEval
tractData_black <- extractdata(modelTract_black, data = data, na.rm = TRUE)
resTract_black <- polr(modelTract_black, data = tractData_black, 
                       method = "probit", Hess = TRUE)
summary(resTract_black)

modelTract_asian <- as.factor(asianHood) ~ asianTract*res_length + 
  edu + income + age + female +latino + black+ democrat + own_home + pocketbookEval
tractData_asian <- extractdata(modelTract_asian, data = data, na.rm = TRUE)
resTract_asian <- polr(modelTract_asian, data = tractData_asian, 
                       method = "probit", Hess = TRUE)
summary(resTract_asian)

modelTract_white <- as.factor(whiteHood) ~ whiteTract*res_length + 
  edu + income + age + female +latino + black+ democrat + own_home + pocketbookEval
tractData_white <- extractdata(modelTract_white, data = data, na.rm = TRUE)
resTract_white <- polr(modelTract_white, data = tractData_white, 
                       method = "probit", Hess = TRUE)
summary(resTract_white)

# Table SI 0.3
# Output was changed by hand for aesthetic purposes
stargazer(resTract_latino, resTract_black, resTract_white, resTract_asian,
          title = "\\textbf{Relationship between Census Measured and Perceived Composition}",
          label = "table2",
          #covariate.labels = labs,
          star.cutoffs = c(0.05, 0.01, 0.001), 
          dep.var.labels = c('Latino', 'Black', 'White', 'Asian American'),
          model.numbers          = FALSE,
          #out = "table_2.tex")
          type = "text")
```

# Residential Length 

```{r, message=FALSE, warning=FALSE}
modelTract_latino <- as.factor(latinoHood) ~ latinoTract + res_length + 
  edu + income + age + female + latino + black + democrat + own_home + pocketbookEval
tractData_latino <- extractdata(modelTract_latino, data = data, na.rm = TRUE)
resTract_latino <- polr(modelTract_latino, data = tractData_latino, 
                        method = "probit", Hess = TRUE)
summary(resTract_latino)

modelTract_black <- as.factor(blackHood) ~ blackTract + res_length + 
  edu + income + age + female +latino + black + democrat + own_home + pocketbookEval
tractData_black <- extractdata(modelTract_black, data = data, na.rm = TRUE)
resTract_black <- polr(modelTract_black, data = tractData_black,
                       method = "probit", Hess = TRUE)
summary(resTract_black)

modelTract_asian <- as.factor(asianHood) ~ asianTract + res_length + 
  edu + income + age + female +latino + black+ democrat + own_home + pocketbookEval
tractData_asian <- extractdata(modelTract_asian, data = data, na.rm = TRUE)
resTract_asian <- polr(modelTract_asian, data = tractData_asian,
                       method = "probit", Hess = TRUE)
summary(resTract_asian)

modelTract_white <- as.factor(whiteHood) ~ whiteTract + res_length + 
  edu + income + age + female +latino + black+ democrat + own_home + pocketbookEval
tractData_white <- extractdata(modelTract_white, data = data, na.rm = TRUE)
resTract_white <- polr(modelTract_white, data = tractData_white, 
                       method = "probit", Hess = TRUE)
summary(resTract_white)

#Table SI 0.2
stargazer(resTract_asian,resTract_black,resTract_latino,resTract_white)
labs <- c("\\% Latino",
          "\\%Black",  
          "\\% White",
          "\\% Asian American",
          "Residential Length",
          "Education","Income","Age","Female",
          "Latino","Black","Democrat","Home Owner",
          "Pocketbook Evaluation")

stargazer(resTract_latino, resTract_black, resTract_white, resTract_asian,
          title = "\\textbf{Relationship between Census Measured and Perceived Composition}",
          label = "table3",
          covariate.labels = labs,
          star.cutoffs = c(0.05, 0.01, 0.001), 
          dep.var.labels = c('Latino', 'Black', 'White', 'Asian American'),
          model.numbers          = FALSE,
          out = "table_3.tex")
```



